Observed2M_Model
Ben Weinstein
June 6, 2016
## [1] "Run Completed at 2016-08-30 23:40:27"#Load in data from Observed2m_Generate.Rmd
load("Observed.RData")
#load("ObservedModel.RData")
newModel<-T#source functions
source("Bayesian/BayesFunctions.R")View Raw Data
ggplot(indat,aes(x=Traitmatch,y=Yobs,col=as.factor(R))) + facet_wrap(~Hummingbird,ncol=4,scales="free") + geom_point() + geom_smooth(method = "glm",method.args=list(family="poisson")) + scale_color_manual("Resource Availability",values=c("black",'blue','red'))
ggplot(indat[indat$Yobs==1,],aes(x=Traitmatch,fill=as.factor(R))) + facet_wrap(~Hummingbird,ncol=4,scales="free") + geom_density(alpha=.7) + scale_fill_manual("Resource Availability",values=c("black",'blue','red'))
m<-glmer(data=indat[,],Yobs~Traitmatch*R+(1|Hummingbird),family="poisson")
summary(m)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Yobs ~ Traitmatch * R + (1 | Hummingbird)
## Data: indat[, ]
##
## AIC BIC logLik deviance df.resid
## 7146.0 7191.3 -3566.0 7132.0 4756
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4850 -0.5293 -0.3050 -0.1925 28.1772
##
## Random effects:
## Groups Name Variance Std.Dev.
## Hummingbird (Intercept) 1.656 1.287
## Number of obs: 4763, groups: Hummingbird, 19
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.079336 0.306095 -3.526 0.000422 ***
## Traitmatch -0.033845 0.004302 -7.867 3.62e-15 ***
## RMedium -0.261394 0.094822 -2.757 0.005839 **
## RHigh -0.303452 0.087808 -3.456 0.000549 ***
## Traitmatch:RMedium -0.002927 0.006711 -0.436 0.662739
## Traitmatch:RHigh 0.001511 0.006168 0.245 0.806440
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trtmtc RMedim RHigh Trt:RM
## Traitmatch -0.154
## RMedium -0.147 0.485
## RHigh -0.161 0.533 0.509
## Trtmtch:RMd 0.101 -0.643 -0.734 -0.347
## Trtmtch:RHg 0.114 -0.711 -0.350 -0.735 0.462indat$pred<-predict(m,type="response")
ggplot(indat,aes(x=Traitmatch)) + geom_point(aes(y=Yobs)) + geom_line(aes(y=pred,col=as.factor(R)))
Hierarchical Model
For hummingbird species i feeding on plant species j observed at time k and sampling event observed by transect (YTransect) or camera (YCamera)
Observation Model:
\[ Yobs_{i,j,k,d} \sim Binomial(N_{i,j,k},\phi) \]
Process Model:
\[ N_{i,j,k} \sim Poisson(\lambda_{i,j,k}) \] \[ log(\lambda_{i,j,k}) = \alpha_i + \beta_{1,i} * Traitmatch_{i,j} + \beta_{2,i} *Resources_k + \beta_{3,i} * Traitmatch_{i,j} * Resources_k \]
Priors
\[ \phi \sim Normal(0,0.386) \] \[\alpha_i \sim Normal(\alpha_\mu,\alpha_\tau)\] \[\beta_{1,i} \sim Normal(\mu_{\beta_1,\tau_{beta_1}})\] \[\beta_{2,i} \sim Normal(\mu_{\beta_2,\tau_{beta_2}})\] \[\beta_{3,i} \sim Normal(\mu_{\beta_3,\tau_{beta_3}})\]
Group Level Logit Transformed Means \[ \mu_\alpha \sim Normal(0,1.67)\] \[\mu_{\beta_1} \sim Normal(0,1.67)\] \[\mu_{\beta_2} \sim Normal(0,1.67)\] \[\mu_{\beta_3} \sim Normal(0,1.67)\]
Group Level Variance \[\tau_{\alpha} \sim Half cauchy(0,1,1)\] \[\tau_{\beta_1} \sim Half cauchy(0,1,1)\] \[\tau_{\beta_2} \sim Half cauchy(0,1,1)\] \[\tau_{\beta_3} \sim Half cauchy(0,1,1)\]
#Source model
source("Bayesian/NmixturePoissonRagged2m.R")
#print model
writeLines(readLines("Bayesian/NmixturePoissonRagged2m.R"))##
## sink("Bayesian/NmixturePoissonRagged2m.jags")
##
## cat("
## model {
## #Compute true state for each pair of birds and plants
## for (i in 1:Birds){
## for (j in 1:Plants){
## for (k in 1:Times){
##
## #Process Model
## log(rho[i,j,k])<-alpha[i] + beta1[i] * Traitmatch[i,j] + beta2[i] * resources[i,j,k] + beta3[i] * resources[i,j,k] * Traitmatch[i,j]
##
## #True State
## S[i,j,k] ~ dpois(rho[i,j,k])
## }
## }
## }
##
## #Observation Model
## for (x in 1:Nobs){
##
## #Observation Process for cameras
## Yobs_camera[x] ~ dbin(detect[Bird[x]],S[Bird[x],Plant[x],Time[x]])
##
## # #Assess Model Fit - Posterior predictive check
## #
## # #Fit discrepancy statistics
## #Camera
## eval_cam[x]<-detect[Bird[x]]*S[Bird[x],Plant[x],Time[x]]
## E[x]<-pow((Yobs_camera[x]-eval_cam[x]),2)/(eval_cam[x]+0.5)
## #
## ynew_cam[x]~dbin(detect[Bird[x]], S[Bird[x],Plant[x],Time[x]])
## E.new[x]<-pow((ynew_cam[x]-eval_cam[x]),2)/(eval_cam[x]+0.5)
## }
##
## #Priors
## #Observation model
## #Detect priors, logit transformed - Following lunn 2012 p85
##
## for(x in 1:Birds){
## #For Cameras
## logit(detect[x])<-dcam[x]
## dcam[x]~dnorm(omega_mu,omega_tau)
## }
##
##
## #Process Model
## #Species level priors
## for (i in 1:Birds){
##
## #Intercept
## #logit prior, then transform for plotting
## alpha[i] ~ dnorm(alpha_mu,alpha_tau)
##
## #Traits slope
## beta1[i] ~ dnorm(beta1_mu,beta1_tau)
##
## #Plant slope
## beta2[i] ~ dnorm(beta2_mu,beta2_tau)
##
## #Interaction slope
## beta3[i] ~ dnorm(beta3_mu,beta3_tau)
##
## }
##
## #OBSERVATION PRIOR
## omega_mu ~ dnorm(0,0.386)
## omega_tau ~ dunif(0,10)
##
## #Group process priors
##
## #Intercept
## alpha_mu ~ dnorm(0,0.386)
## alpha_tau ~ dt(0,1,1)I(0,)
## alpha_sigma<-pow(1/alpha_tau,0.5)
##
## #Trait
## beta1_mu~dnorm(0,0.386)
## beta1_tau ~ dt(0,1,1)I(0,)
## beta1_sigma<-pow(1/beta1_tau,0.5)
##
## #Resources
## beta2_mu~dnorm(0,0.386)
## beta2_tau ~ dt(0,1,1)I(0,)
## beta2_sigma<-pow(1/beta2_tau,0.5)
##
## #Interaction
## beta3_mu~dnorm(0,0.386)
## beta3_tau ~ dt(0,1,1)I(0,)
## beta3_sigma<-pow(1/beta3_tau,0.5)
##
## #derived posterior check
## fit<-sum(E[]) #Discrepancy for the observed data
## fitnew<-sum(E.new[])
##
## }
## ",fill=TRUE)
##
## sink()#Inits
InitStage <- function(){
#A blank Y matrix - all present
initY<-array(dim=c(Birds,Plants,Times),22)
initB<-rep(0.5,Birds)
list(S=initY,dcam=initB)}
#Parameters to track
ParsStage <- c("alpha","beta1","beta2","beta3","alpha_mu","alpha_sigma","beta1_sigma","beta1_mu","beta2_mu","beta2_sigma","beta3_mu","beta3_sigma","detect","E","fit","fitnew")
#Jags Data
Dat<-list(
Yobs_camera = indat$Camera,
Birds=max(indat$jBird),
Bird=indat$jBird,
Plant=indat$jPlant,
Time=indat$jID,
Plants=max(indat$jPlant),
Times=max(indat$jID),
resources=resourceMatrix,
Nobs=nrow(indat),
Traitmatch=jTraitmatch)
#MCMC options
if(newModel){
system.time(
m2<-jags.parallel(data=Dat,parameters.to.save =ParsStage,inits=InitStage,model.file="Bayesian/NmixturePoissonRagged2m.jags",n.thin=10,n.iter=300000,n.burnin=295000,n.chains=2,DIC=F)
)
}## user system elapsed
## 9.944 0.359 14286.424#recompile if needed
runs<-200000
recompile(m2)
if(!newModel){
system.time(m2<-update(m2,n.iter=runs,n.burnin=runs-5000,n.thin=10))
}#extract par to data.frame
pars_detect<-extract_par(m2,data=indat,Bird="jBird",Plant="jPlant")Assess Convergence
###Chains
ggplot(pars_detect[pars_detect$par %in% c("alpha","beta1","beta2","beta3"),],aes(x=Draw,y=estimate,col=as.factor(Chain))) + geom_line() + facet_grid(par~species,scale="free") + theme_bw() + labs(col="Chain") + ggtitle("Species Level")
ggplot(pars_detect[pars_detect$par %in% c("detect"),],aes(x=Draw,y=estimate,col=as.factor(Chain))) + geom_line() + facet_wrap(~species,scale="free",ncol=3) + theme_bw() + labs(col="Chain") + ggtitle("Species Level")
ggplot(pars_detect[pars_detect$par %in% c("beta1_mu","beta1_sigma","beta2_mu","beta2_sigma","beta3_mu","beta3_sigma","alpha_mu","alpha_sigma"),],aes(x=Draw,y=estimate,col=as.factor(Chain))) + geom_line() + theme_bw() + labs(col="Chain") + ggtitle("Group Level Parameters") + facet_wrap(~par,scales="free")
Posteriors
###Posterior Distributions
ggplot(pars_detect[pars_detect$par %in% c("alpha","beta1","beta2","beta3"),],aes(x=estimate)) + geom_histogram(position='identity') + facet_grid(species~par,scales="free") + theme_bw() + ggtitle("Species Parameters")
#Detection figure
pars_detect<-merge(pars_detect,jagsIndexBird,by.x="species",by.y="jBird",all=T)
ggplot(pars_detect[pars_detect$par %in% "detect",],aes(x=par,y=estimate)) + geom_violin(fill='black') + theme_bw() + ggtitle("Detection Probability") + facet_wrap(~Hummingbird)
Detection Table
detecttable<-group_by(pars_detect,Hummingbird,par) %>% filter(par %in% c('detect')) %>% summarize(mean=mean(estimate),lower=quantile(estimate,0.05),upper=quantile(estimate,0.95))
detecttable## Source: local data frame [19 x 5]
## Groups: Hummingbird [?]
##
## Hummingbird par mean lower upper
## (fctr) (fctr) (dbl) (dbl) (dbl)
## 1 Andean Emerald detect 0.2460116 0.13384758 0.3994516
## 2 Booted Racket-tail detect 0.2204931 0.13802587 0.3246381
## 3 Brown Inca detect 0.1666781 0.10561342 0.2510915
## 4 Buff-tailed Coronet detect 0.1877761 0.10501788 0.2843066
## 5 Collared Inca detect 0.1970451 0.10464681 0.3110083
## 6 Crowned Woodnymph detect 0.2204191 0.11881498 0.3418982
## 7 Fawn-breasted Brilliant detect 0.1928675 0.08611181 0.3428303
## 8 Gorgeted Sunangel detect 0.4460895 0.27255985 0.6511495
## 9 Green-crowned Brilliant detect 0.2229545 0.08600007 0.4159489
## 10 Green-fronted Lancebill detect 0.2741602 0.13855144 0.4100603
## 11 Hoary Puffleg detect 0.2360042 0.11421216 0.3988843
## 12 Purple-bibbed Whitetip detect 0.2396576 0.10982839 0.4048972
## 13 Rufous-tailed Hummingbird detect 0.2073415 0.08644013 0.3740578
## 14 Speckled Hummingbird detect 0.1882605 0.10919659 0.2823925
## 15 Stripe-throated Hermit detect 0.2368006 0.13394028 0.3496642
## 16 Tawny-bellied Hermit detect 0.3198216 0.19842256 0.4544295
## 17 Violet-tailed Sylph detect 0.2240017 0.11694238 0.3533654
## 18 Wedge-billed Hummingbird detect 0.1568365 0.07167055 0.2771828
## 19 White-whiskered Hermit detect 0.2468784 0.15660755 0.3671323Discrepancy
The goodness of fit is a measured as chi-squared. The expected value for each day is the detection rate * the estimate intensity of interactions. The expected value is compared to the observed value of the actual data. In addition, a replicate dataset is generated from the posterior predicted intensity. Better fitting models will have lower discrepancy values and be Better fitting models are smaller values and closer to the 1:1 line. A perfect model would be 0 discrepancy. This is unrealsitic given the stochasticity in the sampling processes. Rather, its better to focus on relative discrepancy. In addition, a model with 0 discrepancy would likely be seriously overfit and have little to no predictive power.
fitstat<-pars_detect[pars_detect$par %in% c("fit","fitnew"),]
fitstat<-dcast(fitstat,Draw+Chain~par,value.var="estimate")
ymin<-round(min(fitstat$fit))
ymax<-round(max(fitstat$fit))
ab<-data.frame(x=0:ymax,y=0:ymax)
disc_obs<-ggplot(fitstat,aes(x=fit,y=fitnew)) + geom_point() + theme_bw() + labs(x="Discrepancy of observed data",y="Discrepancy of replicated data",col="Model") + ggtitle("Posterior predictive check") + geom_line(data=ab,aes(x=x,y=y)) + coord_fixed() + ylim(ymin=0,ymax=max(max(c(fitstat$fit,fitstat$fitnew)))) + xlim(xmin=0,xmax=max(max(c(fitstat$fit,fitstat$fitnew))))
disc_obs
#Bayesian p-value
sum(fitstat$fitnew>fitstat$fit)/nrow(fitstat)## [1] 0ggsave("Figures/ObservedDiscrepancy.jpeg",width = 5,height=10)Predicted Relationship
#Expand out pars
castdf<-dcast(pars_detect[pars_detect$par %in% c("beta1_mu","beta2_mu","beta3_mu","alpha_mu"),], Chain + Draw~par,value.var="estimate")Posterior prediction
#Trajectories from posterior
predy<-trajF(alpha=castdf$alpha_mu,beta1=castdf$beta1_mu,trait=indat$Traitmatch,resources=indat$scaledR,beta2=castdf$beta2_mu,beta3=castdf$beta3_mu)
ggplot(data=predy,aes(x=trait)) + geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.4,fill="red") + theme_bw() + ylab("Interactions") + xlab("Difference between Bill and Corolla Length") + geom_point(data=indat,aes(x=Traitmatch,y=Camera)) + geom_line(aes(y=mean)) + geom_point(data=indat,aes(x=Traitmatch,y=Transect)) 
At High and Low Resource Availability
#Trajectories from posterior
predH<-trajF(alpha=castdf$alpha_mu,beta1=castdf$beta1_mu,trait=indat[indat$R=="High","Traitmatch"],resources=indat[indat$R=="High","scaledR"],beta2=castdf$beta2_mu,beta3=castdf$beta3_mu)
predM<-trajF(alpha=castdf$alpha_mu,beta1=castdf$beta1_mu,trait=indat[indat$R=="Medium","Traitmatch"],resources=indat[indat$R=="Medium","scaledR"],beta2=castdf$beta2_mu,beta3=castdf$beta3_mu)
predL<-trajF(alpha=castdf$alpha_mu,beta1=castdf$beta1_mu,trait=indat[indat$R=="Low","Traitmatch"],resources=indat[indat$R=="Low","scaledR"],beta2=castdf$beta2_mu,beta3=castdf$beta3_mu)
predhl<-melt(list(High=predH,Medium=predM,Low=predL),id.vars=colnames(predH))
colnames(predhl)[6]<-"BFlowerL"
predhl$BFlowerL<-factor(as.character(predhl$BFlowerL))
levels(predhl$BFlowerL)<-c("Low","Medium","High")
ggplot(data=predhl,aes(x=trait)) + geom_ribbon(aes(ymin=lower,ymax=upper,fill=BFlowerL),alpha=0.5) + geom_line(aes(y=mean,col=BFlowerL),size=.8) + theme_bw() + ylab("Interactions") + xlab("Difference between Bill and Corolla Length") + geom_point(data=mindat,aes(x=Traitmatch,y=value))+ labs(fill="Resource Availability",col="Resource Availability") 
ggsave("Figures/AllRegression.jpeg",height=5,width=7)Species Predictions
castdf<-dcast(pars_detect[pars_detect$par %in% c("beta1","beta2","beta3","alpha"),], species +Chain +Draw ~par ,value.var="estimate")
#Turn to
castdf$species<-factor(castdf$species,levels=1:max(as.numeric(castdf$species)))
species.split<-split(castdf,list(castdf$species),drop = T)
species.traj<-list()
for(d in 1:length(species.split)){
x<-species.split[[d]]
index<-unique(x$species)
#get data for those species
billd<-indat[indat$jBird %in% index,]
#scale resources
species.traj[[d]]<-trajF(alpha=x$alpha,beta1=x$beta1,beta2=x$beta2,beta3=x$beta3,resources=billd$scaledR,trait=billd$Traitmatch)
}
names(species.traj)<-names(species.split)
species.traj<-melt(species.traj,id.var=colnames(species.traj[[1]]))
#split out names and model
species.traj[,c("Index")]<-colsplit(species.traj$L1,"\\.",c("Index"))
spe<-merge(species.traj,jagsIndexBird,by.x="Index",by.y="jBird")
#plot and compare to original data
ggplot(data=spe,aes(x=trait)) + geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.2,fill='red') + geom_line(aes(y=mean),size=.5) + theme_bw() + ylab("Daily Interaction Rate")+ xlab("Difference between Bill and Corolla Length") + facet_wrap(~Hummingbird,scales="free",ncol=4) + geom_point(data=mindat,aes(x=Traitmatch,y=value),size=2.5) 
Species Predictions: High and Low Availability
castdf<-dcast(pars_detect[pars_detect$par %in% c("beta1","beta2","beta3","alpha"),], species +Chain +Draw ~par ,value.var="estimate")
#Turn to
castdf$species<-factor(castdf$species,levels=1:max(as.numeric(castdf$species)))
species.split<-split(castdf,list(castdf$species),drop = T)
species.traj<-list()
for(d in 1:length(species.split)){
x<-species.split[[d]]
index<-unique(x$species)
#get data for those species
billd<-indat[indat$jBird %in% index,]
sh<-trajF(alpha=x$alpha,beta1=x$beta1,beta2=x$beta2,beta3=x$beta3,resources=billd[billd$R=="High","scaledR"],trait=billd[billd$R=="High","Traitmatch"])
sm<-trajF(alpha=x$alpha,beta1=x$beta1,beta2=x$beta2,beta3=x$beta3,resources=billd[billd$R=="Medium","scaledR"],trait=billd[billd$R=="Medium","Traitmatch"])
sl<-trajF(alpha=x$alpha,beta1=x$beta1,beta2=x$beta2,beta3=x$beta3,resources=billd[billd$R=="Low","scaledR"],trait=billd[billd$R=="Low","Traitmatch"])
sm<-melt(list(High=sh,Medium=sm,Low=sl),id.vars=colnames(sl))
species.traj[[d]]<-sm
}
names(species.traj)<-names(species.split)
species.traj<-melt(species.traj,id.var=colnames(species.traj[[1]]))
#split out names and model
species.traj[,c("Index")]<-colsplit(species.traj$L1,"\\.",c("Index"))
spe<-merge(species.traj,jagsIndexBird,by.x="Index",by.y="jBird")
#plot and compare to original data
indat$rplot<-as.factor(indat$scaledR)
levels(indat$rplot)<-c("Low","Medium","High")
spe$resources<-as.factor(spe$resources)
levels(spe$resources)<-c("Low","Medium","High")
ggplot(data=spe,aes(x=trait)) + geom_ribbon(aes(ymin=lower,ymax=upper,fill=resources),alpha=0.5) + geom_line(aes(y=mean,group=resources),size=.5) + theme_bw() + ylab("Daily Interaction Rate")+ xlab("Difference between Bill and Corolla Length") + facet_wrap(~Hummingbird,scales="free",ncol=4) + geom_point(data=indat,aes(x=Traitmatch,y=Yobs),size=1.2) + scale_fill_manual("Resource Availability",values=c("Grey70","Grey40","Black"))
ggsave("Figures/SpeciesRegression.jpeg",height=7,width=11)Species Level Interaction
castdf<-dcast(pars_detect[pars_detect$par %in% c("beta1","beta2","beta3","alpha"),], species +Chain + Draw~par,value.var="estimate")
#Turn to
castdf$species<-factor(castdf$species,levels=1:max(as.numeric(castdf$species)))
species.split<-split(castdf,list(castdf$species),drop = T)
species.traj<-list()
for(d in 1:length(species.split)){
dat<-species.split[[d]]
index<-unique(dat$species)
#get data for those species
billd<-mindat[mindat$jBird %in% index,]
#Calculate interaction effect
species.traj[[d]]<-intF(alpha=dat$alpha,beta1=dat$beta1,x=billd[billd$value > 0 & !is.na(billd$value),'Traitmatch'],resources=billd[billd$value > 0 & !is.na(billd$value),'scaledR'],beta2=dat$beta2,beta3=dat$beta3)
}
names(species.traj)<-names(species.split)
species.traj<-melt(species.traj,id.var=colnames(species.traj[[1]]))
#split out names and model
species.traj[,c("Index")]<-colsplit(species.traj$L1,"\\.",c("Index"))
spe<-merge(species.traj,jagsIndexBird,by.x="Index",by.y="jBird")
#match colnames
#plot and compare to original data
ggplot(data=spe,aes(x=x)) + geom_ribbon(aes(ymin=lower,ymax=upper,fill=Hummingbird),alpha=0.3) + geom_line(aes(y=mean,col=Hummingbird),size=1) + theme_bw() + xlab("Difference between Bill and Corolla Length") + ylab("Effect of Resources on Trait Difference") + facet_wrap(~Hummingbird,scales="free",ncol=3)
ggsave("Figures/SpeciesInteraction.jpeg",height=6,width=7)These plots can be tricky to interpret if one forgets that trait matching as a covariate is a distance. Therefore remember that a positive slope in the plot above indiciates, “As resources increase species use flowers less similiar to their bill lengths”.
Which species did we predict well?
smat<-pars_detect %>% filter(par=="E") %>% group_by(species) %>% summarize(fit=sum(estimate)) %>% arrange(desc(fit))
smat<-merge(smat,jagsIndexBird,by.x="species",by.y="jBird")Discrepancy matrix
dmat<-pars_detect %>% filter(par=="E") %>% group_by(species,plant) %>% summarize(fit=sum(estimate)) %>% arrange(desc(fit))
dmat<-merge(dmat,jagsIndexBird,by.x="species",by.y="jBird")
dmat<-merge(dmat,jagsIndexPlants,by.x="plant",by.y="jPlant")
head(dmat,10)## plant species fit Hummingbird
## 1 1 4 28.50690 Buff-tailed Coronet
## 2 1 9 35.66932 Green-crowned Brilliant
## 3 1 17 2157.04299 Violet-tailed Sylph
## 4 1 7 921.59691 Fawn-breasted Brilliant
## 5 1 12 30.71657 Purple-bibbed Whitetip
## 6 1 6 60.13242 Crowned Woodnymph
## 7 1 3 8471.06069 Brown Inca
## 8 1 5 663.77643 Collared Inca
## 9 1 10 235.85366 Green-fronted Lancebill
## 10 1 2 69.23956 Booted Racket-tail
## Iplant_Double
## 1 Alloplectus tetragonoides
## 2 Alloplectus tetragonoides
## 3 Alloplectus tetragonoides
## 4 Alloplectus tetragonoides
## 5 Alloplectus tetragonoides
## 6 Alloplectus tetragonoides
## 7 Alloplectus tetragonoides
## 8 Alloplectus tetragonoides
## 9 Alloplectus tetragonoides
## 10 Alloplectus tetragonoides#get order
hord<-dmat %>% group_by(Hummingbird) %>% summarize(n=sum(fit)) %>% arrange(desc(n)) %>% .$Hummingbird
pord<-dmat %>% group_by(Iplant_Double) %>% summarize(n=sum(fit)) %>% arrange(n) %>% .$Iplant_Double
dmat$Hummingbird<-factor(dmat$Hummingbird,levels=hord)
dmat$Iplant_Double<-factor(dmat$Iplant_Double,levels=pord)
ggplot(dmat,aes(x=Hummingbird,Iplant_Double,fill=fit)) + geom_tile() + theme_bw() + scale_fill_continuous("Discrepancy",low="blue",high="red") + theme(axis.text.x=element_text(angle = -90))
Resource Abundance Functions
Let’s take a closer look at distribution of interaction effect posteriors values for each species.
post<-pars_detect %>% filter(par %in% "beta2") %>% group_by(Hummingbird) %>% summarize(mean=mean(estimate),median=median(estimate),lower=quantile(probs=0.025,estimate),upper=quantile(probs=0.975,estimate)) %>% melt(id.vars='Hummingbird')
ggplot(pars_detect[pars_detect$par %in% "beta2",],aes(x=estimate)) + geom_histogram() + facet_wrap(~Hummingbird,scales='free',ncol=4) + geom_vline(data=post,aes(xintercept=value,col=variable))
Interaction density functions
Let’s take a closer look at distribution of interaction effect posteriors values for each species.
post<-pars_detect %>% filter(par %in% "beta3") %>% group_by(Hummingbird) %>% summarize(mean=mean(estimate),median=median(estimate),lower=quantile(probs=0.025,estimate),upper=quantile(probs=0.975,estimate)) %>% melt(id.vars='Hummingbird')
ggplot(pars_detect[pars_detect$par %in% "beta3",],aes(x=estimate)) + geom_histogram() + facet_wrap(~Hummingbird,scales='free',ncol=4) + geom_vline(data=post,aes(xintercept=value,col=variable))
Trait-matching and Bill Length
Do species with long bill lengths have positive traitmatching effects?
#species names
b<-pars_detect[pars_detect$par %in% "beta1",]
#traits
b<-merge(b,hum.morph,by.x="Hummingbird",by.y="English")
post<-b %>% filter(par %in% "beta1") %>% group_by(Hummingbird) %>% summarize(mean=mean(estimate),median=median(estimate),lower=quantile(probs=0.025,estimate),upper=quantile(probs=0.975,estimate),quantile_l=quantile(estimate)[[1]],quantile_u=quantile(estimate)[[2]]) %>% melt(id.vars='Hummingbird')
#get order of mean posterior
ord<-post %>% filter(variable=="mean") %>% arrange(value) %>% .$Hummingbird
b$Hummingbird<-factor(b$Hummingbird,levels=ord)
ggplot(b,aes(y=estimate,x=as.factor(Total_Culmen),group=Hummingbird)) + geom_violin(fill='grey50') + coord_flip() + ggtitle("Trait-matching and Bill Length") + theme_bw()
Interaction and Bill Length
Do species with long bill lengths have positive interaction effects?
#species names
b<-pars_detect[pars_detect$par %in% "beta3",]
#traits
b<-merge(b,hum.morph,by.x="Hummingbird",by.y="English")
post<-b %>% filter(par %in% "beta3") %>% group_by(Hummingbird) %>% summarize(mean=mean(estimate),median=median(estimate),lower=quantile(probs=0.025,estimate),upper=quantile(probs=0.975,estimate),quantile_l=quantile(estimate)[[1]],quantile_u=quantile(estimate)[[2]]) %>% melt(id.vars='Hummingbird')
#get order of mean posterior
ord<-post %>% filter(variable=="mean") %>% arrange(value) %>% .$Hummingbird
b$Hummingbird<-factor(b$Hummingbird,levels=ord)
ggplot(b,aes(y=estimate,x=Hummingbird,fill=Total_Culmen)) + geom_violin() + coord_flip() + scale_fill_continuous(low='blue',high='red') + ggtitle("Interaction Effect and Bill Length") + theme_bw()Estimated niche breadth
castdf<-dcast(pars_detect[pars_detect$par %in% c("beta1","beta2","beta3","alpha"),], species +Chain + Draw~par,value.var="estimate")
#Turn to
castdf$species<-factor(castdf$species,levels=1:max(as.numeric(castdf$species)))
species.split<-split(castdf,list(castdf$species),drop = T)
species.traj<-lapply(species.split,function(dat){
index<-unique(dat$species)
#get data for those species
billd<-indat[indat$jBird %in% index,]
d<-data.frame(alpha=dat$alpha,beta1=dat$beta1,beta2=dat$beta2,beta3=dat$beta3)
#fit regression for each input estimate
sampletraj<-list()
for (y in 1:nrow(d)){
v=exp(d$alpha[y] + d$beta1[y] * billd$Traitmatch + d$beta2[y] * billd$scaledR + d$beta3[y] * billd$Traitmatch*billd$scaledR)
sampletraj[[y]]<-data.frame(x=as.numeric(billd$Traitmatch),y=as.numeric(v),r=as.numeric(billd$scaledR),jBird=billd$jBird,jPlant=billd$jPlant,jTime=billd$jTime)
}
sample_all<-rbind_all(sampletraj)
})
species.traj<-rbind_all(species.traj)Mean Estimates for Corolla Sizes
species.mean<-species.traj %>% group_by(jBird,jPlant,r) %>% summarize(Traitmatch=unique(x),phi=mean(y))
tomerge<-indat %>% select(jBird,jPlant,Hummingbird,Iplant_Double) %>% distinct()
species.mean<-merge(species.mean,tomerge)
#get corolla sizes
species.mean<-merge(species.mean,fl.morph,by.x="Iplant_Double", by.y="Group.1")
#bill order
ord<-hum.morph %>% arrange(Total_Culmen) %>% .$English
species.mean$Hummingbird<-factor(species.mean$Hummingbird,levels=ord)
#add level to hum.morph to match naming convention
species.mean<-merge(species.mean,hum.morph[,c("English","Total_Culmen")],by.x="Hummingbird",by.y="English")
p<-ggplot(species.mean,aes(x=Corolla,y=phi,col=as.factor(r))) + geom_line(size=.9) + geom_vline(aes(xintercept=Total_Culmen),linetype='dashed') + facet_wrap(~Hummingbird,ncol=3,scales="free_y") + theme_bw() + ylab("Probability of Interaction") + scale_color_manual("Resource Availability",labels=c("Low","Medium","High"),values=c("Black","Blue","Red")) + xlab("Flower Corolla Length (mm)")
p
ggsave("Figures/ResponseCurves.jpeg",height=6.5,width=9)Niche Breadth
species.mean<-species.traj %>% group_by(jBird,jPlant,r) %>% summarize(Traitmatch=unique(x),phi=mean(y),phi_low=quantile(y,0.05),phi_high=quantile(y,0.95))
#merge names
species.mean<-merge(species.mean,jagsIndexBird)
species.mean<-merge(species.mean,jagsIndexPlants)
#get corolla sizes
species.mean<-merge(species.mean,fl.morph,by.x="Iplant_Double", by.y="Group.1")
#bill order
ord<-hum.morph %>% arrange(Total_Culmen) %>% .$English
species.mean$Hummingbird<-factor(species.mean$Hummingbird,levels=ord)
#add level to hum.morph to match naming convention
species.mean<-merge(species.mean,hum.morph[,c("English","Total_Culmen")],by.x="Hummingbird",by.y="English")
#label factor
species.mean$r<-as.factor(species.mean$r)
levels(species.mean$r)<-c("Low","Medium","High")
ggplot(species.mean) + geom_ribbon(alpha=0.7,aes(x=Corolla,ymin=phi_low,ymax=phi_high,fill=r)) + theme_bw() + facet_wrap(~Hummingbird,scales="free",ncol=4)+ ggtitle("Niche Breadth") + geom_vline(aes(xintercept=Total_Culmen),linetype='dashed') + geom_line(aes(x=Corolla,y=phi,fill=r)) + ylab("Probability of Interaction") + xlab("Corolla Length (mm)") + scale_fill_manual("Resource Availability",values=c("Grey70","Grey20","Black"))
ggsave("Figures/NicheBreadth.jpeg",height=7,width=11)Range plots
nsplit<-split(species.mean,species.mean$r)
makeR<-function(x){
#input matrix
aggm<-matrix(nrow=nrow(jagsIndexBird),ncol=nrow(jagsIndexPlants),data=0)
for (j in 1:nrow(x)){
aggm[x[j,"jBird"],x[j,"jPlant"]]<-rpois(1,lambda=x[j,"phi"])
}
aggm<-melt(aggm)
colnames(aggm)<-c("jBird","jPlant","P")
tomerge<-species.mean %>% select(jBird,jPlant,Corolla,Hummingbird,Traitmatch) %>% distinct()
aggm<-merge(aggm,tomerge)
return(aggm)
}
nstat<-lapply(nsplit,function(x){
netstat<-melt(lapply(1:100,function(k) makeR(x)),id.vars=c("jBird","jPlant","Hummingbird","Traitmatch","Corolla","P"))
colnames(netstat)[7]<-"Iteration"
return(netstat)
})
names(nstat)<-c("Low","Medium","High")
nstat<-melt(nstat,colnames(nstat[[1]]))Predicted density
#order the levels
nstat$L1<-factor(nstat$L1,levels=c("High","Medium","Low"))
#add level to hum.morph to match naming convention
nstat<-merge(nstat,hum.morph[,c("English","Total_Culmen")],by.x="Hummingbird",by.y="English")
ggplot(nstat[nstat$P>0,],aes(x=Corolla,fill=L1)) + geom_density(alpha=0.8) + facet_wrap(~Hummingbird,scales='free',ncol=4) + theme_bw()+ scale_fill_manual("Resource Availability",values=c("Black","Grey40","Grey85")) + geom_vline(aes(xintercept=Total_Culmen),linetype='dashed')
ggsave("Figures/PredictedDensity.jpeg",dpi=400,height=8,width=10)rangestat<-nstat %>% filter(P==1) %>% group_by(Hummingbird,L1) %>% summarize(mean=mean(Corolla),lower=quantile(Corolla,0.05),upper=quantile(Corolla,0.95))
ggplot(rangestat,aes(x=Hummingbird,ymin=lower,ymax=upper,ymean=mean,col=L1)) + geom_linerange(alpha=0.5,position=position_dodge(width=.75),size=2) + geom_point(aes(y=mean),position=position_dodge(width=.75)) + labs(col="Resource Availability",y="Corolla Length (mm)") + theme_bw() + coord_flip() + theme(legend.position="bottom")
Predicted Mean
meanstat<-nstat %>% filter(P==1) %>% group_by(Hummingbird,L1,Iteration) %>% summarize(a=mean(Corolla))%>% summarize(mean=mean(a),lower=quantile(a,0.05),upper=quantile(a,0.95))
ggplot(meanstat,aes(x=Hummingbird,ymin=lower,ymax=upper,ymean=mean,col=L1)) + geom_linerange(alpha=0.5,position=position_dodge(width=.6),size=2) + geom_point(aes(y=mean),position=position_dodge(width=.6)) + labs(col="Resource Availability",y="Corolla Length (mm)") + theme_bw() + coord_flip() + theme(legend.position="bottom")
Predicted Min
meanstat<-nstat %>% filter(P==1) %>% group_by(Hummingbird,L1,Iteration) %>% summarize(a=min(Corolla))%>% summarize(mean=mean(a),lower=quantile(a,0.05),upper=quantile(a,0.95))
ggplot(meanstat,aes(x=Hummingbird,ymin=lower,ymax=upper,ymean=mean,col=L1)) + geom_linerange(alpha=0.5,position=position_dodge(width=.75),size=2) + geom_point(aes(y=mean),position=position_dodge(width=.75)) + labs(col="Resource Availability",y="Corolla Length (mm)") + theme_bw() + coord_flip() + theme(legend.position="bottom")
Predicted Max
meanstat<-nstat %>% filter(P==1) %>% group_by(Hummingbird,L1,Iteration) %>% summarize(a=min(Corolla))%>% summarize(mean=mean(a),lower=quantile(a,0.05),upper=quantile(a,0.95))
ggplot(meanstat,aes(x=Hummingbird,ymin=lower,ymax=upper,ymean=mean,col=L1)) + geom_linerange(alpha=0.5,position=position_dodge(width=.6),size=2) + geom_point(aes(y=mean),position=position_dodge(width=.6)) + labs(col="Resource Availability",y="Corolla Length (mm)") + theme_bw() + coord_flip() + theme(legend.position="bottom")
Generate network
#Split by resource
nsplit<-split(species.mean,species.mean$r)
makeN<-function(x){
#input matrix
aggm<-matrix(nrow=nrow(jagsIndexBird),ncol=nrow(jagsIndexPlants),data=0)
for (j in 1:nrow(x)){
aggm[x[j,"jBird"],x[j,"jPlant"]]<-rpois(1,lambda=x[j,"phi"])
}
#calculate network statistic
nstat<-networklevel(aggm,index=c("niche overlap"),level="lower")
}
nstat<-lapply(nsplit,function(x){
netstat<-melt(t(sapply(1:500,function(k) makeN(x))))
colnames(netstat)<-c("Iteration","Metric","value")
return(netstat)
})
names(nstat)<-c("Low","Medium","High")
nstat<-melt(nstat,colnames(nstat[[1]]))
nstat$L1<-factor(nstat$L1,c("Low","Medium","High"))
levels(nstat$Metric)<-"Niche Overlap"
ggplot(nstat,aes(x=value,fill=L1)) + geom_density(alpha=0.7) + facet_wrap(~Metric,scales='free',nrow=3) + scale_fill_discrete("Resource Availability") + theme_bw() + scale_fill_manual("Resource Availability",values=c("Grey80","Grey50","Black"))
ggsave("Figures/NetworkStatistics.jpeg",height=4,width=6,dpi=600) Compared to raw visits
sindat<-split(indat,indat$scaledR)
ndat<-lapply(sindat,function(x){
web<-acast(x[x$Yobs==1,],Hummingbird~Iplant_Double,value.var = "Yobs", fun= function(x){ length(unique(x))})
nstat<-t(networklevel(web,index=c("niche overlap"),level="lower"))
})
names(ndat)<-c("Low","Medium","High")
ndat<-melt(ndat)
colnames(ndat)<-c("Var1","Metric","value","L1")
ggplot(nstat,aes(x=value,fill=L1)) + geom_density(alpha=0.6) + facet_wrap(~Metric,scales='free',nrow=2) + scale_fill_discrete("Resource Availability") + geom_vline(data=ndat,aes(xintercept=value,col=L1),linetype="dashed") + scale_color_discrete(guide='none') 
save.image("ObservedModel.RData")